Versuch einer Objektdefinition über Minima¶
Vielleicht ist es möglich Objekte aus dem IR-10.8-µm-Kanal über die Minima und davon ausgehenden Grenzen zu definieren.
Um das zu untersuchen, brauchen wir zunächst ein paar Pakete.
In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage as ndi
from skimage.feature import peak_local_max
import cv2
import datetime as dt
import l15_msevi.msevi as msv
import sys
sys.path.append("/vols/talos/home/stephan/utils/tracking")
import tracking_common as tc
Beispieldaten laden¶
Wir laden für den 28.07.2012, 12:00 UTC ein Beispiel.
In [2]:
t = dt.datetime(2012,7,28,12,0)
s = msv.MSevi(time=t,chan_list=["IR_108"])
s.lonlat()
s.rad2bt()
ir108 = s.bt['IR_108']
Nun suchen wir nach den lokalen Minima im Bild.
In [18]:
minima = peak_local_max(-ir108, min_distance=2,threshold_abs=-273)
In [9]:
print minima
In [19]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
x = [m[1] for m in minima]
y = [m[0] for m in minima]
ax.scatter(x,y,marker='+',color='blue')
Out[19]:
Das sind ziemlich viele Minima. Es gibt mehrere Möglichkeiten deren Anzahl zu reduzieren:
- das Eingangsbild glätten und maskieren
- den Mindestabstand vergrößeren
- den Schwellwert strenger wählen.
Das probieren wir im Folgenden mal aus.
Glättung und Maskierung des Eingangsbildes¶
Für die Glättung gibt es mehrere Möglichkeiten. Hier probieren wir Folgende aus:
- Mittelwertglättung
- Gauß'scher Filter
- Medianglättung
- bilateraler Filter
All diese Verfahren sind abhängig von der Filterbreite. Was ist ein guter Wert dafür?
In [4]:
ir108_transformed = (tc.scale_array_min_max(ir108)*255).astype("uint8")
mean_blur = cv2.blur(ir108_transformed,(5,5))
gauss_blur = cv2.GaussianBlur(ir108_transformed,(5,5),0)
median_blur = cv2.medianBlur(ir108_transformed,5)
bilateral_blur = cv2.bilateralFilter(ir108_transformed,15,15,15)
In [49]:
plt.imshow(ir108_transformed)
plt.colorbar()
Out[49]:
Darstellung der Glättungsergebnisse¶
In [53]:
fig,ax = plt.subplots(2,2,figsize=(14,14))
ax[0,0].imshow(mean_blur,vmin=0,vmax=255,cmap='gray_r')
ax[0,0].set_title(u"Mittelwertglättung")
ax[0,1].imshow(gauss_blur,vmin=0,vmax=255,cmap='gray_r')
ax[0,1].set_title(u"Gaußglättung")
ax[1,0].imshow(median_blur,vmin=0,vmax=255,cmap='gray_r')
ax[1,0].set_title(u"Medianglättung")
ax[1,1].imshow(bilateral_blur,vmin=0,vmax=255,cmap='gray_r')
ax[1,1].set_title(u"bilateraler Filter")
Out[53]:
Die Glättung mit dem bilateralen Filter sieht subjektiv am besten aus. Das Feld ist relativ stark geglättet, aber die Kanten sind noch gut erhalten.
In [115]:
thresh = (273 - np.min(ir108)) / (np.max(ir108) - np.min(ir108))*255
In [116]:
thresh
Out[116]:
In [117]:
minima = peak_local_max(-bilateral_blur, min_distance=2,threshold_abs=thresh)
In [118]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
x = [m[1] for m in minima]
y = [m[0] for m in minima]
ax.scatter(x,y,marker='+',color='blue')
Out[118]:
Das sieht schon besser aus. Nur die Häufungen von Minima sind seltsam.
Vergrößerung des Mindestabstandes¶
In [111]:
minima = peak_local_max(-ir108, min_distance=5,threshold_abs=-273)
In [112]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
x = [m[1] for m in minima]
y = [m[0] for m in minima]
ax.scatter(x,y,marker='+',color='blue')
Out[112]:
Das sieht besser aus.
Strengere Schwellwerte¶
In [113]:
minima = peak_local_max(-ir108, min_distance=2,threshold_abs=-250)
In [114]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
x = [m[1] for m in minima]
y = [m[0] for m in minima]
ax.scatter(x,y,marker='+',color='blue')
Out[114]:
Das sind immernoch viele Minima, aber sinnvoller verteilt.
Kombination aus allem¶
In [5]:
th = 250
thresh = (th - np.min(ir108)) / (np.max(ir108) - np.min(ir108))*255
test = 255 - bilateral_blur
thr = 255 - thresh
minima = peak_local_max(test, min_distance=15,threshold_abs=thr)
#from skimage.morphology import extrema
#minima = extrema.h_minima(bilateral_blur,30)
In [6]:
minima
Out[6]:
In [7]:
test = 255 - bilateral_blur
In [12]:
from mpl_toolkits.axes_grid1 import make_axes_locatable
def colourbar(mappable):
ax = mappable.axes
fig = ax.figure
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
return fig.colorbar(mappable, cax=cax)
In [13]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
#ax.imshow(np.ma.masked_less(minima,1))
test_plot = ax.imshow(test,vmin=0,vmax=255,cmap='gray')
colourbar(test_plot)
Out[13]:
In [9]:
print(np.max(ir108) - np.min(ir108))
print(th - np.min(ir108))
print thresh
In [14]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
x = [m[1] for m in minima]
y = [m[0] for m in minima]
ax.scatter(x,y,marker='+',color='blue')
Out[14]:
In [15]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
test_plot = ax.imshow(np.ma.masked_less(test,thr))
colourbar(test_plot)
Out[15]:
In [190]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
test_plot = ax.imshow(np.ma.masked_greater(ir108,250))
colourbar(test_plot)
Out[190]:
Das sieht schon recht gut aus. Bloß dies Häufung von Minima ist seltsam. Das kann man wahrschenlich folgendermaßen beseitigen:
- Matrix mit Nullen erzeugen
- Punkte die ein Minimum haben auf 1 setzen
- diese Matrix markieren
- Massenschwerpunkte der verschiedenen Objekte bestimmen und als neue Minima betrachten
In [16]:
from scipy import ndimage as ndi
new_minima = np.zeros_like(ir108)
for m in minima:
new_minima[m[0],m[1]] = 1
lab,nlab = ndi.measurements.label(new_minima)
new_mins = ndi.center_of_mass(new_minima,lab,np.arange(1,nlab))
In [17]:
new_mins
Out[17]:
In [19]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ir_plot = ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
colourbar(ir_plot)
x = [m[1] for m in new_mins]
y = [m[0] for m in new_mins]
ax.scatter(x,y,marker='+',color='blue')
Out[19]:
Das sieht als repräsentativer für die Objekte aund Startpunkte für die Objektdefinition ziemlich gut aus. Es sind Minima von nichtkonvektiven Wolken dabei, aber um diese zu entfernen, ist wahrscheinlich eine Wolkentypisierung nötig.
In [27]:
ct = msv.MSGtools.get_nwcsaf_prod('CT',t,calibrate=True)
In [28]:
ct = np.ma.masked_where(np.logical_or(ct<5, ct>10),ct)
In [29]:
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ct_plot = ax.imshow(ct)
colourbar(ct_plot)
Out[29]:
